import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
import sympy as sp
from sympy import symbols, lambdify, integrate, diff

# 设置字体
plt.rcParams['font.family'] = 'DejaVu Sans'
plt.rcParams['axes.unicode_minus'] = False
plt.rcParams['font.size'] = 10

def lagrange_interpolation(x_data, y_data):
    """
    Lagrange Interpolation Method
    """
    n = len(x_data)
    x = symbols('x')
    
    # 构造拉格朗日基函数
    def l_i(i):
        product = 1
        for j in range(n):
            if j != i:
                product *= (x - x_data[j]) / (x_data[i] - x_data[j])
        return product
    
    # 构造拉格朗日插值多项式
    P = sum(y_data[i] * l_i(i) for i in range(n))
    
    return P, x

def newton_interpolation(x_data, y_data):
    """
    Newton Interpolation Method
    """
    n = len(x_data)
    x = symbols('x')
    
    # 计算差商
    def divided_differences():
        # 初始化差商表
        dd = [y_data.copy()]
        
        for i in range(1, n):
            row = []
            for j in range(n - i):
                numerator = dd[i-1][j+1] - dd[i-1][j]
                denominator = x_data[j+i] - x_data[j]
                row.append(numerator / denominator)
            dd.append(row)
        
        return dd
    
    # 计算差商
    dd = divided_differences()
    
    # 构造牛顿插值多项式
    P = dd[0][0]  # f[x0]
    
    for i in range(1, n):
        term = dd[i][0]  # f[x0, x1, ..., xi]
        for j in range(i):
            term *= (x - x_data[j])
        P += term
    
    return P, x

def best_square_approximation():
    """
    Best Square Approximation for f(x) = |x| in S = span{1, x², x⁴}
    """
    x = symbols('x')
    
    # 定义基函数
    phi0 = 1
    phi1 = x**2
    phi2 = x**4
    
    # 目标函数 f(x) = |x|
    f = abs(x)
    
    # 计算内积 <f, phi_i> 和 <phi_i, phi_j>
    def inner_product(f1, f2):
        return integrate(f1 * f2, (x, -1, 1))
    
    # 构建法方程组
    A = np.zeros((3, 3))
    b = np.zeros(3)
    
    # 计算矩阵元素
    A[0,0] = float(inner_product(phi0, phi0))  # <1,1> = 2
    A[0,1] = float(inner_product(phi0, phi1))  # <1,x²> = 2/3
    A[0,2] = float(inner_product(phi0, phi2))  # <1,x⁴> = 2/5
    A[1,0] = A[0,1]  # 对称矩阵
    A[1,1] = float(inner_product(phi1, phi1))  # <x²,x²> = 2/5
    A[1,2] = float(inner_product(phi1, phi2))  # <x²,x⁴> = 2/7
    A[2,0] = A[0,2]  # 对称矩阵
    A[2,1] = A[1,2]  # 对称矩阵
    A[2,2] = float(inner_product(phi2, phi2))  # <x⁴,x⁴> = 2/9
    
    # 计算右端项
    b[0] = float(inner_product(f, phi0))  # <|x|,1> = 1
    b[1] = float(inner_product(f, phi1))  # <|x|,x²> = 1/2
    b[2] = float(inner_product(f, phi2))  # <|x|,x⁴> = 1/3
    
    # 求解线性方程组
    coeffs = np.linalg.solve(A, b)
    
    # 构造逼近多项式
    a, b_coeff, c = coeffs
    p = a + b_coeff * x**2 + c * x**4
    
    return p, x, coeffs

def solve_problem_1():
    """
    Problem 1: Lagrange and Newton Interpolation
    """
    print("=" * 60)
    print("Problem 1: Lagrange and Newton Interpolation")
    print("=" * 60)
    
    # 给定数据点
    x_data = [1, -1, 2]
    y_data = [0, -3, 4]
    
    print(f"Data points: {list(zip(x_data, y_data))}")
    
    # 拉格朗日插值
    print("\n1. Lagrange Interpolation:")
    P_lagrange, x = lagrange_interpolation(x_data, y_data)
    print(f"Lagrange polynomial: {P_lagrange}")
    print(f"Simplified: {sp.simplify(P_lagrange)}")
    
    # 牛顿插值
    print("\n2. Newton Interpolation:")
    P_newton, x = newton_interpolation(x_data, y_data)
    print(f"Newton polynomial: {P_newton}")
    print(f"Simplified: {sp.simplify(P_newton)}")
    
    # 验证结果
    print("\n3. Verification:")
    P_lagrange_func = lambdify(x, P_lagrange, 'numpy')
    P_newton_func = lambdify(x, P_newton, 'numpy')
    
    for i, (xi, yi) in enumerate(zip(x_data, y_data)):
        lagrange_val = P_lagrange_func(xi)
        newton_val = P_newton_func(xi)
        print(f"Point {i+1}: x={xi}, y={yi}")
        print(f"  Lagrange: P({xi}) = {lagrange_val}")
        print(f"  Newton: P({xi}) = {newton_val}")
        print(f"  Difference: {abs(lagrange_val - newton_val):.2e}")
    
    # 绘制图像
    x_plot = np.linspace(-2, 3, 100)
    y_lagrange = P_lagrange_func(x_plot)
    y_newton = P_newton_func(x_plot)
    
    plt.figure(figsize=(12, 8))
    
    plt.subplot(2, 2, 1)
    plt.plot(x_plot, y_lagrange, 'b-', linewidth=2, label='Lagrange Interpolation')
    plt.plot(x_data, y_data, 'ro', markersize=8, label='Data Points')
    plt.xlabel('x')
    plt.ylabel('y')
    plt.title('Lagrange Interpolation Result')
    plt.grid(True, alpha=0.3)
    plt.legend()
    
    plt.subplot(2, 2, 2)
    plt.plot(x_plot, y_newton, 'r-', linewidth=2, label='Newton Interpolation')
    plt.plot(x_data, y_data, 'ro', markersize=8, label='Data Points')
    plt.xlabel('x')
    plt.ylabel('y')
    plt.title('Newton Interpolation Result')
    plt.grid(True, alpha=0.3)
    plt.legend()
    
    plt.subplot(2, 2, 3)
    plt.plot(x_plot, y_lagrange, 'b-', linewidth=2, label='Lagrange', alpha=0.7)
    plt.plot(x_plot, y_newton, 'r--', linewidth=2, label='Newton', alpha=0.7)
    plt.plot(x_data, y_data, 'ro', markersize=8, label='Data Points')
    plt.xlabel('x')
    plt.ylabel('y')
    plt.title('Comparison of Both Methods')
    plt.grid(True, alpha=0.3)
    plt.legend()
    
    plt.subplot(2, 2, 4)
    diff_plot = y_lagrange - y_newton
    plt.plot(x_plot, diff_plot, 'g-', linewidth=2)
    plt.xlabel('x')
    plt.ylabel('Difference')
    plt.title('Lagrange - Newton Difference')
    plt.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('problem1_interpolation_english.png', dpi=300, bbox_inches='tight')
    plt.show()
    
    return P_lagrange, P_newton

def solve_problem_2():
    """
    Problem 2: Best Square Approximation
    """
    print("\n" + "=" * 60)
    print("Problem 2: Best Square Approximation")
    print("=" * 60)
    
    # 求解最佳平方逼近
    p, x, coeffs = best_square_approximation()
    
    print(f"Best approximation polynomial: p(x) = {p}")
    print(f"Coefficients: a = {coeffs[0]:.6f}, b = {coeffs[1]:.6f}, c = {coeffs[2]:.6f}")
    
    # 计算逼近误差
    f = abs(x)
    error_squared = (f - p)**2
    error_integral = integrate(error_squared, (x, -1, 1))
    error_norm = float(sp.sqrt(error_integral))
    
    print(f"L2 norm error: {error_norm:.6f}")
    
    # 绘制对比图
    x_plot = np.linspace(-1, 1, 1000)
    f_plot = np.abs(x_plot)
    p_func = lambdify(x, p, 'numpy')
    p_plot = p_func(x_plot)
    
    plt.figure(figsize=(15, 10))
    
    plt.subplot(2, 2, 1)
    plt.plot(x_plot, f_plot, 'b-', linewidth=2, label='f(x) = |x|')
    plt.plot(x_plot, p_plot, 'r--', linewidth=2, label='Best approximation p(x)')
    plt.xlabel('x')
    plt.ylabel('y')
    plt.title('Best Square Approximation Result')
    plt.grid(True, alpha=0.3)
    plt.legend()
    
    plt.subplot(2, 2, 2)
    error_plot = f_plot - p_plot
    plt.plot(x_plot, error_plot, 'g-', linewidth=2)
    plt.xlabel('x')
    plt.ylabel('Error')
    plt.title('Approximation Error f(x) - p(x)')
    plt.grid(True, alpha=0.3)
    
    plt.subplot(2, 2, 3)
    plt.plot(x_plot, f_plot**2, 'b-', linewidth=2, label='f^2(x)')
    plt.plot(x_plot, p_plot**2, 'r--', linewidth=2, label='p^2(x)')
    plt.xlabel('x')
    plt.ylabel('y^2')
    plt.title('Squared Functions Comparison')
    plt.grid(True, alpha=0.3)
    plt.legend()
    
    plt.subplot(2, 2, 4)
    squared_error = (f_plot - p_plot)**2
    plt.plot(x_plot, squared_error, 'm-', linewidth=2)
    plt.xlabel('x')
    plt.ylabel('Squared Error')
    plt.title('Squared Error (f(x) - p(x))^2')
    plt.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('problem2_approximation_english.png', dpi=300, bbox_inches='tight')
    plt.show()
    
    # 计算一些特定点的值
    print(f"\nSpecific point values:")
    test_points = [-1, -0.5, 0, 0.5, 1]
    for xi in test_points:
        f_val = abs(xi)
        p_val = p_func(xi)
        error = abs(f_val - p_val)
        print(f"x = {xi:4.1f}: f(x) = {f_val:6.3f}, p(x) = {p_val:6.3f}, Error = {error:6.3f}")
    
    return p, coeffs, error_norm

def main():
    """
    Main function
    """
    print("Numerical Analysis Problems Solution")
    print("=" * 60)
    
    # 解决第一题
    P_lagrange, P_newton = solve_problem_1()
    
    # 解决第二题
    p, coeffs, error_norm = solve_problem_2()
    
    print("\n" + "=" * 60)
    print("Summary")
    print("=" * 60)
    print("1. Lagrange and Newton interpolation give identical results")
    print("2. Best square approximation solved successfully with small error")
    print("3. All plots saved as PNG files")
    
    return P_lagrange, P_newton, p, coeffs

if __name__ == "__main__":
    main()
